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The effects of finite size and of finite current on the current- voltage characteristics of Josephson 
junction arrays are studied both theoretically and by numerical simulations. The cross-over from 
non-linear to linear behavior at low temperature is shown to be a finite size effect and the non-linear 
behavior at higher temperature, T > Tkt, is shown to be a finite current effect. These are argued 
to result from competition between the length scales characterizing the system. The importance of 
boundary effects is discussed and it is shown that these may dominate the behavior in small arrays. 
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Two dimensional (2D) Josephson junction arrays 
(J J A) in the absence of an external magnetic field have 
been extensively studied Jlj-Q as model systems to ob- 
serve the Kosterlitz-Thouless (KT) transition §,§. In 
the thermodynamic limit when the array size N — > oo 
and the applied current I — > such that IN — > oo the 
current-voltage (/V") relation is V ~ 7 a ( T ) where the ex- 
ponent a(T) > 3 when the temperature T < Tkt and 
a(T) = 1 when T > T KT @0k However, it is observed 
in experiments on small arrays - !]!]] and in simulations fj^] 
that the IV relation for T < Tkt becomes linear when 
the applied current / is lowered below some threshold 
value. An explanation that this is due to a residual mag- 
netic field inducing some free vortices was proposed ffl] 
and an opposing point of view that this is an intrinsic 
effect was offered [12[. However, experiments on large ar- 
rays g show no cross-over to linear behavior at T < Tkt 
down to the lowest accessible currents which implies that 
a low T linear IV relation is a finite size effect. In this ar- 
ticle, we demonstrate both analytically and numerically 
that this is a finite size effect and that the nonlinear IV 
relation becomes linear for / < I x « I /N where / is the 
applied current per junction and I is the critical cur- 
rent. It is also observed that the IV relation is also non- 
linear at temperatures slightly above Tkt and we show 
that this is a finite current effect for / > If w I /£+(T) 
where £+(T) is the dimensionless correlation length for 
T > T KT - 

The behavior of the system is controlled by the length 
scales in the problem: the linear size N of the array, the 
current length £j = I D /I which is the maximum size of 
a bound vortex pair in the presence of of a current / 
and the thermal correlation length £(T) which diverges 
as T — > T^ T and is infinite for T < Tkt- Above Tkt, 
£ = £ + (T) may be interpreted as the size of the largest 
bound pair which is stable against thermal fluctuations, 
while bound pairs of all sizes exist for T < Tkt- The 
essential difference is that, for T > Tkt, there is also 
a finite density of thermally excited free vortices lead- 
ing to a linear resistivity R ~ £+ 2 while for T < Tkt 
the only free vortices are due to the current unbinding 
of vortex pairs which leads to a non-linear IV relation. 
In the 2DXY model which describes the system, there is 



a fourth length scale £-(T) which also diverges at Tkt 
and £ + ~ £ 27r This length plays a rather different 

role to the others as it is a measure of the scale at which 
the renormalized vortex interaction K(l) is essentially at 
its asymptotic value K(oo) = Kr(T). A renormalization 
group analysis shows that £_ is the scale at which devi- 
ations from critical behavior become significant and not 
the scale at which the RG equations become invalid. The 
other scales N, are scales at which the RG equa- 

tions do become invalid and, to proceed, some physically 
motivated approximation such as Debye-Hiickel [JHJ must 
be made and the result extrapolated back. 

A scaling assumption for the resistivity R which is con- 
sistent with the renormalization group is 



V/I = e-' l K(Ne- l ,€e- l ,Zre- l ,l/ln4-) 



(1) 



where the dynamical exponent z = 2 and 1Z is an un- 
known scaling function. This scaling form is familiar 
from more conventional scaling functions with the excep- 
tion of the combination I / which is a consequence of 
the marginally relevant and irrelevant scaling fields in the 
system. The ultimate effect of this is to produce temper- 
ature and current dependent power laws. One can obtain 
some information from the scaling ansatz by choosing a 
value of I at which the scaling function 1Z can be com- 
puted. For example, in the case < N, £, one can choose 
e l = at which scale all vortex pairs are unbound by 
the current and the vortices may be regarded as moving 
independently in a viscous medium driven by the applied 
current so that 

V/I = I 2 n(N/^,^/^, (2) 

For T < Tkt, N — > oo and In^i > ln£^ we expect 

K(oo, oo, 1, Inti/lnt-) - /^-(T) _ jx(t) (3) 

where £_(T) = exp(l/x(T)) with x(T) = b\T-T KT \ 1/2 , 
which leads to the standard nonlinear IV relation 

V ~ jkK r (T)+i ag ^ e renorma ii zec i stiffness constant 
-kKr(T) = 2 + x(T) and, when ln£i < we expect 

V ~ I 3 . Finite size dominated behavior is also contained 
in Eq.(0) by taking N < £/ and e l — N to obtain 
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V/I = N- 2 K{\, oo, ii/N, InN/ln^i) (4) 

Now, when InN > ln£_, identical arguments lead to a 
linear IV relation V/I ~ n-*Kr(t) and yyj _ 
when ZniV < ln£— . By making appropriate assumptions 
about the behavior of the unknown scaling function TZ in 
the various limits, all the expected behaviors can be re- 
produced and one expects qualitative changes when one 
of the ratios of lengths is of order unity. Of course, real 
dynamical calculations should be done to determine the 
functional form of TZ at the scale / and then use the RG 
equations to extrapolate to physical values of the param- 
eters. We have been unable to do this explicitly but hope 
to address this in the future. The strategy of this paper is 
to numerically simulate the IV relation and to interpret 
the data by these scaling considerations. 

We have performed simulations on unfrustrated square 
N x N arrays of various sizes from N — 4 to N — 64 at 
different temperatures to measure the voltage across the 
array as a function of the external current / per bond. 
We use a modification of the Langevin dynamical method 
HQ as proposed by Falo, Bishop and Lomdahl Q . We 
assume that all junctions have the same critical current 
I and are shunted by equal resistances R Each super- 
conducting grain at the vertices of the lattice has a ca- 
pacitance C to ground. The dynamical equations for the 
phase 9 n and the voltage V n of the superconducting grain 
at site n follow from charge conservation and the Joseph- 
son equation. 

d9 n /dt = 2eV n /h 
CdV n /dt = 4 E sin(9 m - 9 n ) + 

<m> 

R' 1 E(^™-^)+ E *%n (5) 

<m> <m> 

The sum over < m > is over the nearest neighbors of 
site n and I^ n is a thermal noise current in the bond mn 
satisfying the fluctuation dissipation relation 

< I%{t)I%{t') >= (2T/R)(S ik S jl - S tl S 3k )S(t - t') (6) 

The capacitances C to ground are introduced purely for 
computational convenience and the simulations were per- 
formed with an intermediate value of the McCumber- 
Stewart parameter Q = 2e 2 R 2 I C/h = 1. 

To imitate a typical experimental configuration, we 
use current injection and extraction from superconduct- 
ing busbars at two opposite edges of the array. The 
busbars are modelled by columns of N superconduct- 
ing grains coupled by infinitely strong junctions, each 
of these grains again having a capacitanc C to ground so 
that, in the absence of a magnetic field, the phase 9 and 
the potential V will be the same everywhere on a bar. 
The busbars are each connected to the array by a set of 
N junctions which are assumed identical to those in the 
array. With this geometry, separate equations of motion 
are necessary for the phases 9l,9r and potentials Vl,Vr 



of the busbars on the left and right edges of the array. 
Modelling the bars in this way, they are very similar to 
the equations of motion for the phases and potentials of 
the sites in the interior of the array. 

d9 L /dt = 2eV L /h 
d9 R /dt = 2eV R /h 

N 

CdV L /dt = I + IoN- 1 E sin{9i -6 L ) + 

i=l 

N 

R^N-^iVi-VL) 

i=l 

N 

CdV R /dt = -I + IoN- 1 E sin(0i - Or) + 

i=l 

N 

BT^-^iVi-Vn) (7) 

i=l 

where / is the current per bond injected into the left bus- 
bar and extracted from the right. The sums over i are 
over the N sites connected to the busbars. 

This busbar geometry has some advantages over the 
more usual method of uniform current injection where 
an equal current is injected or extracted from each site 
at the edges of the array and is not technically difficult to 
implement. In fact, it is easier in the presence of an exter- 
nal magnetic field as the gauge invariant phase is constant 
along the busbars and the currents in the bonds attached 
to them are automatically properly accounted for, in con- 
trast to uniform injection. The potential drop across a 
piece of the array in the presence of an applied current is 
due to free vortices, which may be produced by thermal 
excitation, by the applied current or by an external mag- 
netic field, being driven across the array perpendicular to 
the current. Uniform injection emphasizes dissipation at 
the edges where the current is injected and extracted as 
vortices are easily created there and will be attracted to 
the edges by their images. Unbinding of a vortex/image 
pair is not needed for dissipation but vortices at the edges 
are driven across the array by the applied current and, in 
a finite system, most of the dissipation may occur at the 
edges and not in the bulk. This implies that there is a sig- 
nificant voltage drop across narrow regions at the edges 
and that these boundary effects may dominate in finite 
arrays. In the busbar geometry used in the simulations 
reported here, vortices cannot be created adjacent to a 
superconducting busbar and are repelled from the array 
edges so that dissipation there will be suppressed. One 
thus expects that there will be a negligible voltage drop 
at the edges and all will occur in the bulk of the array. 
We have done simulations in both geometries to confirm 
this picture and find that, with the system sizes used, 
most of the potential drop is at the edges with uniform 
injection. 
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FIG. 1. Phase difference <j>(t) across a N = 32 array at 
T = 0.8 and I/I = 0.025. 1 oooc 



Simulations were performed on N x N square ar- 
rays in the busbar geometry with N — 4, 8, 16, 32, 64 
at T = 0.8 . 1. 0, 1.1, 1.3 with the critical temperature > 0.0100 , 
T KT fa 0.89 @ where T is in units of M /2e. The equa- 
tions of motion, Eqs.(^-^|), were integrated using the sim- 
ple Euler algorithm with a time step At = 0.05 in units 00010 r 
of 1/uj = {hC/2eI o y/ 2 , the inverse Josephson plasma 
frequency. Decreasing the time step by a factor of 10 
did not change the results. The system is initially in a 0001 o.oi 
configuration at t = with all 9 n = = V n and inte- 
grated over at least 10 6 time steps and the phase differ- 
ence 4>(t) = 6ji{t) — 0L(t) recorded as a function of t (see , . 000 o 
Fig.l). From this, it is clear that (f>(t) grows by discon- 
tinuous jumps of magnitude 2im due to the motion of 
vortices. The mean voltage drop V across the system is ° 1 °oo r 
given by 



V/RI = m r )-<t>{fi))/tr (8) <: 



0.0001 

0.01 



where t r is the run time in units of l/u)j. To estimate 
the errors, we divide each run into four equal intervals 
and estimate V for each interval to obtain AV/V w 0.1. 
Since the voltage V is caused by vortices crossing the 
array, the number n v of these is 



2im v = <f>(t r ) - 0(0) (9) 



The run time t r is chosen so that n v > 100 and one ex- 0.001 

— 1/2 

pects an error in V of AV/V ~ n v < 0.1 which is 
consistent with the estimate from block averaging. 



FIG. 2. IV relation (i = I/I versus v = V/RI ) 
with current injection from busbars for arrays of sizes 
TV = 4, 8, 16, 32, 64. Errors are about the size of symbols. 
(a)T = 0.8; (b)T = 1.0; (c)T = 1.1; (d)T = 1.3 

The results of the simulations are summarized in 
Figs.2(a-d) where the IV relation at fixed T is shown 
for several system sizes. In Fig. 2(a) the data is shown 
for T = 0.8 < Tkt- At this temperature, we expect that 
£_ is small relative to all array sizes N and to £/ so that 
the voltage V will depend mainly on the ratio £//iV and, 
from our scaling hypothesis of Eq.(|l|), we expect that 

V ~ C7° (T) when £/ < N and V ~ I when & > N with a 
crossover between the two behaviors when £i/N = O(l). 
The data of Fig. 2(a) is completely consistent with this 
hypothesis with the cross-over from non-linear to linear 
Ohmic behavior occurs at a current I x /I = 1/N corre- 
sponding to = N. This cross-over current is in good 
agreement with the experimental data on small arrays: 
see Figs. 3 and 10 of M. The theoretical explanation of 
this cross-over to Ohmic behavior at small currents at 
T < Tkt is quite simple. The origin of the non-linear 
IV relation in 2D superconductors at T < Tkt is that 
vortex pairs of separation r > £/ = I Q /I are unbound 
by the applied current / at a rate ~ (7 / ' j o ^k r {t) anc j 
reformed at rate ~ n^. In the steady state, these rates 
are equal so that n v ~ (I / ' I o y KRtyT \ Since the resistiv- 
ity V/I cx n v , this yields the AHNS result |J] V ~ I a{T) 
with a(T) = nK R (T)+l. We note that the IV relation of 
Fig. 2a agrees fairly well with the AHNS value for the ex- 
ponent a(T) at T = 0.8 when the theoretical estimate is 
a(T) w 3.8. However, our simulations are not consistent 
with a recent proposal that a(T) — 2ttKji(T) - la 4.5 
at T = 0.8 Q|. If £/ > A, no vortex pairs are unbound 
by the applied current as there are no pairs of separation 
r > . so that the only source of free vortices is by ther- 
mal excitation in the finite system with n v ~ N~ vKr ^ 
which yields a small linear resistivity which decreases as 
N increases. 

When T > Tkt, the thermal length scale is £+( r ) ~ 
(^{-t) H where r = (T — T KT )/T K t is the reduced 
temperature [ [i"3"| . This correlation length is much larger 
than £-(|r|) for the same |r| and cannot be considered 
small relative to or to N for all our arrays. At very 
small applied currents I / 1 < l/£+, the system contains 
bound vortex pairs of separation r < and also free 
vortices of areal density n v oc £+ 2 . The dissipation will 
be dominated by these free vortices leading to a linear 
Ohmic IV relation with a resistivity oc n v . However, at 
larger currents £j < £+ the current length scale will con- 
trol the dissipation and one expects a similar mechanism 
of vortex pair unbinding and recapture as when r < 
leading to V - I a{T ^ with 1 < a(T) < 3. These qual- 
itative arguments imply that the crossover regime from 
the large current power law TV relation to the final lin- 
ear Ohmic relation for small I / 1 is very wide, possibly 
several orders of magnitude when T « Tkt and N > 
which is consistent with experimental observation 



The simulation results for T > Tkt are shown in 
Figs. 2(b), (c),(d). From these, it is clear that the cross- 
over from a linear to a non-linear IV rlation is governed 
by competition between three length scales N, £/, £+(T). 
For T = 1.3 (Fig. 2d) where we expect the thermal length 
to be rather small, the IV relation is independent of 
the array size and the data may be interpreted by as- 
suming £ + < N for all arrays and that the cross-over 
from non- linear to linear behavior is controlled by / £j 
with £+(T = 1.3) ~ 4. At the two lower temperatures 
of Figs. 2(b), (c), the data may be similarly interpreted 
with the added complication of additional finite size in- 
duced cross-overs. At T — 1.0, Fig. 2(b) shows that the 
three largest systems with N — 16, 32, 64 start to deviate 
from power law behavior V ~ I a at the same current //, 
which we interpret as indicating that / £/ is the con- 
trolling quantity. The two smaller arrays with N = 4, 8 
deviate from the others when £i/N > 1, roughly as they 
do when T < Tkt- We see both from numerical simula- 
tions and scaling arguments that the linear IV relation 
in superconducting arrays below Tkt for small currents 
is entirely due to finite size effects and the non-linear 
behavior above Tkt is a finite current effect. A con- 
sistent interpretation of the data of Figs. 2(b), (c),(d) for 
T > Tkt may be made and we infer that the correla- 
tion lengths £+(T) ~ 16,8,4 at T = 1.0,1.1,1.3 in fair 
agreement with estimates from equilibrium simulations 
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FIG. 3. IV relation (i = I/I versus v = V/RI ) with 
uniform current injection for sizes N = 4, 16, 64 at T = 0.8. 

All the simulations discussed so far were performed in 
the busbar geometry and with free boundary conditions 
in the transverse direction to imitate the standard ex- 
perimental configuration. However, several simulations 
have been done with uniform injection and extraction of 
the current and with periodic boundary conditions in the 
transverse direction. As discussed earlier, uniform injec- 
tion is expected to emphasize dissipation at the edges of 
the array. To compare with current injection from bus- 
bars, we also did simulations at T — 0.8 with uniform 
injection and transverse periodic BC for three array sizes 
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N = 4, 16,64 which are shown in Fig. 3. Particularly for 
N = 64, the TV relation is linear at much larger current 
than for the busbar geometry which implies that there is 
some extra dissipation. For N — 4, the crossover seems 
to occur at about the same current as in the busbar geom- 
etry. The explanation of this is that vortices are readily 
nucleated at the edges as vortex/image pairs and, once 
created, are attracted to the edge so that much of the 
voltage drop or dissipation occurs not in the bulk but at 
the edges. This is confirmed by studying the phase pro- 
file across the array. In Fig. 4(a) is shown the phase pro- 
file for injection from busbars. The phase is constant in 
thin layers adjacent to the bars and all the phase change 
occurs in the bulk. Since the voltage drop AV oc A9 
from Eq.(j|) there is no dissipation at the edges and it is 
all in the bulk. On the other hand, Fig. 4(b) shows the 
phase profile across the array for uniform injection. As 
expected, the major part of the phase change occurs in 
two thin layers at the boundaries which corresponds to 
motion of free vortices along the edges. This will give 
rise to a linear IV relation with a finite resistance which, 
although an edge effect which presumably is negligible in 
the thermodynamic limit, will dominate in a finite system 
at sufficiently low current. 
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FIG. 4. Phase profile of array at time 
t r . N = 64,/// = 0.05, T = 0.8. The average phase of 
the grains in the n th column plotted against column number 
n. (a)Current injection from busbars at n = 0, 65 with free 
BC in transverse direction; t r — 5.10 5 . (b)Uniform injection; 
t r — 2.10 5 and periodic BC in transverse direction. 

In conclusion, this work shows that experimental and 
simulation studies of the IV characteristics of supercon- 
ducting arrays are liable to be dominated by finite size 
and edge effects unless great care is taken to minimize 
them. Large array sizes seem to be essential to make 
any meaningful comparison between theory and experi- 
ment or numerical simulation. When performing a DC 
measurement of the IV relation, dissipation at the edges 
should be minimized by current injection from busbars 
as edge effects can give a large contribution to the lin- 
ear resistivity which may dominate the non-linear bulk 
contribution. Hopefully, this paper gives some indication 
how to take account of such effects and will motivate re- 
search into finding a complete scaling form for the voltage 
for finite size and finite currents to obtain a better un- 
derstanding of the complete IV curves and not just in 
special limits as we do at present. 
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